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Abstract 

A general orthogonal coordinate system is used to describe various 
axisymmetric and two-dimensional shapes. Close approximations to planetary 
probe configurations are possible. The full Navier-Stokes equations are 
discretized in this coordinate system in a manner based on Allen and 
Cheng's numerical procedure. The bow shock is treated as a discontinuity 
which floats between grid points. Completely coupled flows over the 
forebody, base, and near wake have been calculated over a cylinder, 
sphere, and an approximation to the Viking Aeroshell for 2 10 and 

1000 £ Re £ 30,000. The program gives good comparisons with experimental 
data. A major contribution of this work is that it allows one to study 
the effect of changes in the blunt body shape on the base flow structure. 
Also, some problem areas in determining the base flow for increasing 
Reynolds number are discussed. In particular, it is found that the mean 
free path of the fluid near the wall immediately below the corner of the 
Viking Aeroshell, which experiences a severe expansion, can become greater 
than the local mesh size required to resolve the boundary layer in the fore- 
body. Negative pressures, densities, and temperatures can be calculated in 
these instances. It appears that part of this problem is due to an 
inability to capture a lip shock close to the wall with the given grid. 


*Aero-Space Technologist, Aerothermodynamics Branch, Space Systems 
Division. Presently involved in Graduate Study at Princeton University, 
Princeton, New Jersey. 
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Nomenclature 

transformations constants, nondimensionalized by Rn^ 
speed of sound > m/sec 

coordinate stretching constants used in Eq. (3) 

shock shape parameters in Eq. (6) 

functions defined in Eq. (3) 

total enthalpy nondimensional ized by V ^ 

metric coefficient, nondimensional ized by Rj^ 

indices on node points (Fig. 2) 

heat transfer coefficient, W/m-°K 

length of recirculation region along axis, nondimensional ized 
by R„ 

Mach number 

total number of node points in 0 direction 
total number of node points in n direction 
direction normal to wall, m 
Prandtl number 

O 

pressure nondimensional ized by 
2 

heat transfer, W/m 

Reynolds number - based on nose radius of curvature 
nose radius of curvature, m 
temperature, nondimensional ized by 
freestream temperature, °K 
stagnation temperature, °K 


t 


time nondimensional ized by R../V 
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u,v velocities directed along lines of constant r and 

e respectively, nondimensionalized by 
freestream velocity, m/sec 

x,y,z Cartesian coordinates, nondimensionalized by 

Y metric, nondimensionalized by 

= 1 (2D case) 

= distance from axis in axi symmetric case 

Y ratio of specific heats 

6 shock standoff distance, normal to body nondimensionalized by R| 

initial guess on shock standoff distance on an axis, 
nondimensionalized by Rj^ 
e smoothing constant in Eq. (5) 

n transformed r coordinate 

0,r,(J) transformed coordinates 

0ggp value of 0 where flow separates in the base 

p viscosity, nondimensionalized by y(T^) 

Y Mach angle 

p density nondimensionalized by 

p^ freestream density, kg/m^ 

^ angle defined by Eq. (7) and Fig. 3 

Subscripts : 
w wall 


oo 


freestream 
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Introduction 

Currently projected missions to probe the atmospheres of the outer 
planets have raised questions concerning the influence of base flow on 
the low density aerodynamics of a probe and the radiative heating levels 
imposed on a probe afterbody. The radiative heating levels on the 
forebody will be severe and although the ratio of afterbody to forebody 
heating is expected to be small the magnitude of the afterbody heating 
can still be significant. A lip shock standing off the body just 
upstream of the separation point adds to the complexity of the base flow, 
and makes it difficult to solve the base flow problem on a computer. 

There are two general approaches for attacking this problem 
computationally: 1) The coupling of several approximate techniques 

which are designed to model specific features of the flowfield, the 

so called "patchwork" approach, and 2) the solution of the full Navier- 

12 3 4 

Stokes equations in the entire flowfield. ’ ’ ’ Both of these 
approaches have their own advantages and disadvantages. Several 
approximate techniques for the solution of wake flows have been 
presented in the 1 iterature^’^’^*®’^ and summaries of these results 
are available. Solutions to the Navier-Stokes equations in the 
base region using specified conditions on the inflow boundary have also 
been presented in the literature. In general, these approximate 
techniques offer the advantage of obtaining solutions rapidly on the 
computer. However, the nature of the various approximations can distort 
the results and usually these techniques are restricted to specific body 
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geometries and flow regimes. A solution of the full Navier-Stokes equations 
in the entire flowfield eliminates the need for making the various 
assumptions required for the approximate solutions. However, solutions 
of the Navier-Stokes equations typically require a great deal of 
computational time and storage. Ideally, one should be able to use the 
two solution approaches to complement each other; i.e., use the Navier- 
Stokes solutions to improve the nature of the approximations in the 
approximate techniques and use these rapidly running, specialized 
techniques to perform parametric studies for design analyses. 

The material presented herein is a continuation of a research 
effort^ which is directed toward developing a computer code which can 
compute full Navier-Stokes solutions around complete probe-like 
configurations for a wide range of supersonic Mach numbers and Reynolds 
numbers. As Mach number and Reynolds number are increased various flow 
phenomena develop (shock strengths increase; boundary layers, lip shocks, 
and recompression shocks form) which require special consideration 
within the computational code. Flowfield solutions at low supersonic 
Mach numbers and low Reynolds numbers were described in the previous 
work and good agreement was found with experimental data. Whereas the 
technique presented in Ref. 1 employed shock capturing throughout 
the flowfield, a procedure for treating the bow shock as a discontinuity 
which is allowed to float within the computational grid has been included 
in the present analysis. This capability permits the computation of 
higher Mach number and Reynolds number flows. An improved coordinate 
stretching capability has extended the Reynolds number range by making 



it possible to pack more points in the boundary layer. Shock and body 
slip conditions have been included in the present code. Observations 
of flowfield characteristics and comparisons to experimental data are 
presented for a Mach number range 2 < _< 10 and a Reynolds number 

range 1000 1 Rg 1 30,000. In keeping with a philosophy of adding only . 
one complicating factor at a time, only laminar, perfect gas, 
two-dimensional, or axisymmetric flows are considered. 

It is also noted that this program is still considered in a 
developmental stage. Although comparisons with experimental data are 
encouraging there are still aspects of this approach which should be 
investigated more fully (these potential problem areas will be discussed 
in later sections). The purpose of the present paper is to report the 
progress of this investigation and to point out some interesting 
problem areas and flow phenomena that have been encountered. 
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Coordinate System 

1 5 

A generalized curvilinear orthogonal coordinate system which can 
be used for approximating various axisymmetric and two-dimensional body 
shapes; including spheres, ellipses, spherically-capped cones, flat-faced 
cylinders with rounded corners, circular disks and planetary probes; is 
used herein. The transformation from the (6, r, 4>) domain to the 
(x, y, z) domain for an axisymmetric coordinate system is written as: 


N 


x(6,r,4>) = (-B sinh r + C cosh r) cos 6-2] A e"^ cos n0 

n =2 " 




N 


y(6»r,(l)) = [(B cosh r - C sinh r) sin 0 + n0]cos 4» ) (1) 


n=2 


N 


z(0,r,(j)) = [(B cosh r - C sin r) sin 0 + n0]sin (P 

ri=2 " 




where N is a positive integer greater than two and A^, B, and C 
are arbitrary constants. A two-dimensional transformation to the x,y 
plane is obtained by setting ({) = 0. Lines of constant r are 
transformed to circles in the x,y plane as r increases without limit 
through negative values. The line r = 0 forms the body in the x,y 
plane. The outflow boundary of the computational space is mapped to 
infinity (i.e. a circle of infinite radius) using an additional 
transformation. An example of a body with the associated coordinate 
system is shown in Fig. 1. The metric coefficients h0 and h^ are 
equal in this coordinate system, thus saving considerable computational 
time and storage. 



Analytic Development 


The Navier-Stokes equations written in an orthogonal curvilinear 
coordinate system may be obtained from the literature.^ ^ Since only 
the steady-state solution is desired, a term involving the temporal 
derivative of pressure has been omitted from the energy equation in order 
to simplify the numerical procedure. The convective terms in the 
governing equations are written in conservation form as recommended^^ 
for shock capturing whereas the dissipative terms are expanded in 
terms of u, v, H, and y, (Even though the bow shock is fitted it 
is assumed that we will be able to capture the lip shock and 
recompression shock in the base.) 

In preliminary numerical calculations it was found that finite- 

difference approximations of terms involving products of metric 

coefficients with the conservative flow variables (i.e. 8(hpu)/80, 

2 2 

3[h (p+pu )]/39 etc.) caused large errors in the numerical solution 
at large distances from the body where coordinate stretching was 
significant. Writing the governing equations in a form which separated 
derivatives of metric coefficients, which are evaluated analytically, 
from derivatives of conservation flow variables, which are evaluated 
numerically, eliminated this problem. The governing equations expanded 
in this manner become: 


9p/3t = "[3pu/30 + 3pV/3r + pu(l/Y 3Y/39 + 1/h 3h/30) 


+ pv(l/Y 3Y/3r + 1/h 3h/3r]/h a 

3pu/3t = -[3(p + pu^)/30 + 3puv/3r - pv^ 1/h 3h/30 

+ pu^(l/Y 3Y/30 + 1/h 3h/30) + puv(l/Y 3Y/3r + 2/h 3h/3r)]/h 
+ F^/(h\) b 

3pv/3t = -[3puv/30 + 3(p + pv^)/3r - pu^ 1/h 3h/3r 

+ pv^(l/Y 3Y/3r + 1/h 3h/3r) + puv(l/Y 3Y/30 + 2/h 3h/30)]/h 
+ c 

3pH/3t = -[3puH/30 + 3pvH/3r + puH(l/Y 3Y/30 + 1/h 3h/30) 

+ pvH(l/Y 3Y/3r + 1/h 3h/3r)]/h 

+ p3/(h\) + F4/(h\p^) d 

where F -|_4 are various dissipation functions defined in Ref. 1 with 
y = y(T) using Sutherland's law and Pr = constant. 

An additional coordinate transformation is utilized which simplifies 
the treatment of the outflow boundary by mapping it to infinity where 
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conditions are known and gives some control in the density of mesh points 
near the body. The new coordinate n is defined so that 

a 

b ^(3) 
c 

This stretching does not affect the orthogonality of coordinate system. 
Derivatives with respect to r are rewritten as derivatives with 
respect to ^ as follows; 


r = -Hn{f^)/f 2 

f^ = 1 + a^O - n) 
^3 

f2 = 320 + a^n 


8f/3r = 8f/9n 3ri/9r 


3^f/3r^ = 3^f/3n^(3n/3r)^ + 3f/3n 3^n/3r^ 


where 3n/9r and 


3^n/9r^ 


can be obtained algebraically. 



Numerical Approach 


All of the material presented in this section, except for the 
discussion on shock floating and initial conditions, is described in 
Ref. 1. For the sake of brevity and completeness, only a brief summary of 
material presented in the earlier work is presented. 

A modification of the Brailovskaya scheme which was introduced by 
13 18 

Allen and Cheng ’ was applied to the governing equations. This 
particular method was chosen because the viscous stability limit 
in the numerical procedure is removed. The stability limit on the time 
step for the inviscid terms is written as 


At < l/[|u|/(hA0) + |v|/(h 9r/9n An) 


+ a 




1/(A0^) 


+ l/(9r/9ri Aii 



(4) 


where a is the local speed of sound. 

Preliminary calculations indicated that a smoothing routine was 

necessary to eliminate numerical instabilities which occur in the 

vicinity of a captured shock and in the wake. A nonphysical damping 

1 9 

function was used to eliminate these instabilities. Terms of fourth 
order in the spatial grid are used to smooth the variables p, pu, pv, 
and pH after every iteration. For example 



- e[(A0)^9V80^ + (An)^ 9^p/8n^ 

* sJ 



(5) 



12 


The tilde indicates the undamped results from the corrector step 
of the difference scheme and e is a constant such that 0 < e £ 1/24. 

Small variations to the smoothing formulas must be implemented near the 
boundaries j = 1 , j = NJ and in the vicinity of the fitted bow shock 
(Fig. 2). 

Initial Conditions - The initial shock shape is prescribed by the 
formula (similar to one by Van Dyke^^) 

4 - ( 6 ) 

where 

B, = 1/(1-M^) , 0 < c < 2, R > 0 

x-j = X - Xq, and Xq ~ ^ x(tt,0) 

(The constant c = 1 yields an hyperbola, c 0 causes the Mach angle, 

V, to be approached more quickly, and Rs controls the magnitude of 
yg at a particular x location). The points of intersection (values 
of r) where lines of constant 0 intersect the shock are determined 
implicitly. The shock angle is calculated and conditions on the shock 
at these points are evaluated using the Rankine-Hugoniot equations. The 
no-slip conditions are imposed on the body. A stagnation pressure is 
calculated and the pressure distribution around the body is assumed to 
vary as the cosine squared of the body angle, until the freestream 
pressure is achieved, at which point the pressure on the body is held 
at p^. An isentropic expansion around the body is used to approximate 


the distribution of p and H. The quantities p, p, u, and v are 
assumed to vary linearly in n between the shock and the body along 
lines of constant 6 , and H is calculated from the equation of state. 

A linear variation of these quantities in n based on freestream values 
of p, p, u, and v is used in the base region where lines of 
constant 0 do not intersect the shock. 

The components of the uniform freestream velocity field in the 
present coordinate system are calculated as v = - sin 41 and 
u = cos 4j, where 4j is the angle of the vector tangent to a line of 
constant r in a direction of increasing 6 , measured in a counter- 
clockwise direction from the horizontal (Fig. 3). The differentials 

dx and dy along a line of constant r can be determined from 
Eq. (1) and may be used to show that 


cos \p = 


sin ip = 


N 

(B sinh r - C cosh r) sin 0 + X) nA e*^*^ sin n0 '' 

n =2 ” 

h 

N 

(B cosh r - C sinh r) cos 0 + X) nA cos n 0 f 

n =2 " 


a 

b 


Shock Floating - The shock floating technique, treating the shock as 
discontinuity which lies between mesh points in the computational field 
(Fig. 2), has been described by Moretti Salas, and Daywitt, 

pA 

et al . Treating the bow shock as a discontinuity using the Rankine- 
Hugoniot equations eliminates the problem, associated with shock 
capturing, of creating a numerical overshoot and undershoot of properties 
near the shock, possibly causing static enthalpy to go negative. The 
choice of shock floating over nondimensional izing by shock displacement 
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was made for two reasons. Shock displacement along a line of constant 0 
becomes unbounded in the base region; and, for small values of 6 near 
the base symmetry plane, lines of constant 0 never intersect the shock. 
Also, nondimensional izing by shock displacement destroys the orthogonality 
of the coordinate system and consequently further complicates the 
governing equations. 

The shock position is tracked along lines of constant 0. This 
coordinate is the obvious choice in the stagnation region because lines 
of constant 0 are nearly normal to the bow shock here. However, in 
sweeping around the body in a clockwise direction, it is found that 
these lines intersect the shock at smaller angles; consequently, a small 
displacement of shock position in a direction normal to the shock (which 
is the direction shock velocities are determined) is seen as a large 
displacement along a line of constant 0. The shock position could be 
tracked along lines of constant n in these situations; however, the 
bookkeeping and internal logic would become more complex. It has been 
found that lines of constant 0 can be used to track the shock movement 
so long as one ensures that the shock angle never becomes less than the 
Mach angle during the iteration process. 

The shock floating algorithm cannot be described in detail here 
because of space limitations but a general description is supplied below. 
The numerical procedure at interior points is described in Ref. 1. 

Special difference formulas are used when the shock cuts a computational 
molecule. A predictor step is performed at mesh points behind and on the 
bow shock. Pressure behind the shock and the shock angle determines shock 
velocity. These velocities are resolved along coordinate directions to 
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obtain shock displacement. The Rankine-Hugoniot equations for a moving 
shock are used to redefine p, u, v, and H on the shock. The 
procedure is repeated for the corrector step. Some damping of the shock 
motion and smoothing of the shock shape are necessary early in the 
iteration process; however, these fixes can be turned off when the shock 
approaches a steady position. The shock converges rapidly in the 
stagnation region and more slowly downstream where the shock angle 
approaches the Mach angle. Depending on the severity of the "misalignment" 
of the initial shock shape, a wave in the shock can be observed to travel 
downstream. The use of central differences to calculate shock angle in 
the forebody region, where the flow behind the shock is subsonic or 
transonic, and forward differences downstream, prevented these waves from 
reflecting back into the stagnation region, thus speeding convergence. 
Sometimes this wave caused the shock angle to become less than the Mach 
angle. In such situations, rather than use the Rankine-Hugoniot equations 
which are invalid, shock displacement is defined in a manner to force 
the shock angle to equal the Mach angle. 
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Results and Discussion 

All of the material presented herein, unless otherwise noted, was 
computed on a grid of 51 x 50, with a^ = 1,5, a2 = 1 . = 2, and 

a^ = 1.5. The constants describing the cylinder are B = C = 1, = 0. 

The constants describing the Viking Aeroshell (Fig. 1) are B = 1.5146477, 

C = 0.71607873, A2 = -0.059786238, A3 = -0.056648540, and A4 = 0.042817883. 

General Observations - For the given initial conditions, a region 
of high density and pressure builds up in the base area and eventually 
flows out toward the outflow boundary. The formation of this high density 
region produces gradients which can cause instabilities in the wake. 

These instabilities can be controlled by adjusting the value of the 
smoothing parameter, e, and by decreasing the integration time step to 
some fraction, {» 1/10) of the stability limit of Eq. (4). At times, 
it is necessary to iterate only in the base flow region, while holding 
the rest of the field constant. These problems become more severe as 
and Re^ are increased. Many of these problems can be alleviated 
by starting the procedure with a low Reynolds number for the first few 
thousand iterations. It is advisable to iterate the solution during 
this period in blocks of 500 iterations in order to ascertain the effect 
of changes in the smoothing parameter, the time step, or the Reynolds 
number, always picking up in the new solution where the old solution 
left off. When one has passed the hurdle of integrating in time beyond 
these artificially induced transients and instabilities the time step 
can be raised to the full stability limit and e can be decreased. 
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Note that many of these problems could be overcome by establishing 
more realistic initial conditions in this coordinate system. Also, it 
is possible that a given set of initial conditions may result in 
transients that are too complex to be resolved by the present technique. 
This is especially true for the initial guess on shock shape at high 
Mach numbers. The smoothing that was incorporated into the program in 
the vicinity of the shock is designed to allow the solution to progress 
from a relatively large error in the initial shock shape. Still, it is 
possible to specify an initial shock shape that cannot be iterated to 
convergence. In such cases it is useful to study the history of the 
shock motion (i.e.. Is the shock moving in or out? Is the shock bluntness 
in the nose region increasing or decreasing?) and adjust the parameters 
of the initial shock shape (Eq. (6)) accordingly. It should also be 
noted that the fix of using a low Reynolds number initially, as mentioned 
in the previous paragraph, can sometimes make the shock floating process 
more difficult (i.e., there really is no discrete shock to track). In 
such cases some compromise must be found in choosing the values of 
Re and e to be used. Experience dictates the best choice and is the 
best teacher in these situations. 

Convergence is achieved within 15,000 iterations. It is defined by 
looking at the value E = - p*^)/p*^| and at the distribution 

of properties along the base symmetry line. (When the maximum value of 
E in the field becomes less than lO""^ and properties on the base line 
of symmetry are constant to three decimal places after 1000 iterations 
then the solution is said to be converged.) After a converged solution 
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has been obtained for a particular Mach number and Reynolds number, the 
solution at a new Reynolds number can be obtained within approximately 
5,000 iterations using the old solution as an initial condition. The 
same is true for changes in that produce small changes in v. 

Some distortion of properties near the boundary n = 0 is observed. 
The poor resolution makes it impossible to measure vehicle drag based 
on an integrated momentum deficit in the wake. Mapping the outflow 
boundary to infinity has proved useful in these early calculations in 
the sense that no stability problems have ever been encountered due to 
this particular treatment of the boundary. However, it is believed 
that a more exact treatment of the outflow boundary at some finite 
distance downstream will make it possible to obtain better resolution 
of the wake throat. 

Comparisons with Experimental Data : Comparisons between the present 

25 

method and experimental results of Tewfic and Giedt for pressure 
distribution and heat transfer on a cylinder are presented in Figs. 4 
and 5. Shock and body slip conditions were used although the 
calculated effects of shock slip were negligible. Also, it appears the 
no-slip pressure distribution gives a better agreement to experimental 
data than the slip condition in the stagnation region (Fig. 4) though 
the effect is small. A detailed comparison of base pressure 
distributions agree well with experimental data for similar flow 
conditions. The flow separation point in Fig. 4 occurs slightly 
downstream of the pressure minimum, at an inflection point in the curve. 
The heat transfer results for the M =5.5 case (Pr = 0.77 
corresponding to 90° K wall temperature) agree well with the experimental 
data, except in the stagnation region where there is approximately an 
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8-percent difference (Fig. 5). The = 5.73 case yield values of 
q which agree well with experimental data in the stagnation region, 
fall below experimental data for much of the forebody, and level off 
slightly above experimental data in the base. In both cases, the 
calculated heat transfer was obtained according to the formula 

q = k 3T/8n + yu 8u/8n (dimensional equation) 

or 

where the second term is due to "sliding friction" , which accounts for 
slip flow effects. Possible explanations for these discrepancies 
include the assumptions of Pr = constant, the use of Sutherland's 
law at these temperatures (T^ ~ 40*^ K) and experimental error. Changes 
in the heat transfer results were negligible when the smoothing parameter, 
e, was varied from e = 0.003 to e = 0.0005. 

The static pressure distribution along the wake centerline of an 
adiabatic cylinder is compared to experimental data of McCarthy and 
Kubota in Fig. 6. The distributions near peak centerline pressure 
are also plotted for the previous cases to illustrate the effect of 
Reynolds number. The overprediction of pressure downstream is most 
likely due to lack of resolution in the numerical technique. The 
experimental data for a laminar wake show a dip in the transverse 
pressure distribution across the centerline which cannot be resolved 
in the present computational grid. Velocity vectors around the cylinder 
are shown in Fig. 7. Although there is no distinct line where the 
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velocity vectors are sharply turned due to the presence of a recompression 

shock one can see the gradual turning of the vectors to a parallel 

direction in the vicinity of where the shock should be, even in the 

relatively coarse grid. (The distribution of node points can also be 

understood from this figure. Every other node point, in a checkerboard 

pattern, has a velocity vector associated with it, the arrows tail 

originating at the node point.) The discrepancy of base pressure is 

27 

within the experimental accuracy (Ap (e)/p (it) = ±.02). 

W W 

Wake centerline pressure distributions and velocity vectors are 

presented in Figs. 6 and 8 for a sphere at the conditions given above 

for = 5.64. The recirculation region is seen to be slightly 

larger for the sphere. The maximum recirculating centerline velocity 

and Mach number for the sphere is greater than that of the cylinder 

^''max "" v^ax "" 0*^25 for cylinder, = 0.436 for sphere, 

^max ~ 0.288 for cylinder). Excellent agreement for forebody pressure 

distribution and shock shape was obtained between the present method 

28 

and the numerical method of Graves (Fig. 9). 

Viking Aeroshell - Calculations have been made for supersonic 
flow over an approximation to the Viking Aeroshell in C02*» y = 1.285, 

Pr = 0.685, adiabatic wall; for the following conditions: M^ = 2, 

Re^ = 1000, 5000, 30,000; M = 5, Re = 5000; M = 10, Re = 1000, 5000. 

Fig. 10 presents the isobars for the M = 2, Re = 1000 case. The 

oo 00 

arrows and dashed lines indicate the shift in position of the isobars 
for Re^ = 30,000. The expansion fan and the coalescence of isobars 
to form the recompression shock are evident in this figure. Drag 
coefficients and some properties of the near wake flow (separation 


point, Qggp. length of recirculation region, L, base pressure and 

temperature at 6 = o, P,^(o)/p^, Ty^(o)/T^, and maximum pressure on wake 

centerline, P^a^/Poo) presented in Table U 

The recirculation pattern for the = 2, Re^ = 5000 case is 

presented in Fig. 11. There are two recirculating flow patterns 

(clockwise rotation) separated by a small counterclockwise circulating 

flow. The upper recirculating region in this case seems to form a 

horizontal platform over which the separated flow passes. Separation 

occurs slightly downstream of the upper corner of the body corresponding 

29 13 

to observations made by Hama and Allen and Cheng for a sharp corner 
body. Increasing the Mach number while holding Re^ constant causes 
the upper recirculation region to decrease in size and the separation 
flow angle decreases below the horizontal. At = 10 this upper 
recirculation region vanishes. Supersonic recirculation velocities on 
the axis of symmetry were calculated for M^ = 2, Re^ = 1000, 5000, 
30,000 and M = 5, Re = 5000. For M = 2, Re = 1000 sonic velocity 
was just achieved and as the Reynolds number was increased a Mach number 
of approximately 1.4 was attained and a small shock formed approximately 
0.4 nose radii away from the base, with a radius of approximately 0.25 
nose radii. The recirculation shock strengths for M^ = 5, Re^ = 5000, 
and M^ = 2, Re^ = 1000 are approximately equal, suggesting that such 
phenomena are a function of M^/Ae^, though clearly more calculations 
are needed. A calculation by Erdos and Zakkay for supersonic flow 
over a wedge in which nearly supersonic recirculating velocity was 
reported is the only other mention of large recirculating velocity that 
was found in the literature. They conclude that the recirculation 
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velocity varies directly with the vorticity entering the region, up to 
some choking condition in the flow. The present results are consistent 
with that observation. Several numerical experiments in which e and grid 
size were varied did nothing to change the position or strength of the 
shock, although changes in resolution caused some differences in pressure 
downstream (Fig. 12). The calculation of this "recirculation shock" is 
a disturbing development, even though the magnitudes of pressure and 
density are small in its vicinity. It has been assumed that the base 
flow is steady and axisymmetric, and possibly the imposition of this 
symmetry at the axis artificially produces the calculated phenomenon. 
Whatever the cause, whether it is truly physical or numerical, it is 
an interesting result that certainly needs more investigation before the 
question of its existence can be satisfactorily resolved. 

A calculation was attempted for = 5, Re^ = 30,000. This case 
could not be run to convergence because negative values for density, 
pressure, and enthalpy were calculated near the wall just below the 
corner of the probe. The expansion around the corner was extremely 
severe. Even a half cell away from the body the velocity was supersonic. 
Increased resolution near the body did not alter this situation. In 
fact, a normal mesh spacing that was just fine enough to resolve the 
boundary layer in the forebody was less than the mean free path 
immediately behind the expansion corner. No discernible boundary layer 
is present in this region. There is a rapid compression at the wall 
following the expansion, the flow becomes subsonic and then separates. 

This behavior may be due to the formation of a lip shock, normal to the 
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body, immediately below the expansion corner. This shock could not be 
captured with the present distribution of node points in a tangential 
direction and with the present formulation of the boundary conditions. 

The numerical undershoot of properties on the low pressure side of the 
shock that occurs with shock capturing would then explain the calculation 
of negative enthalpies. Indeed, it may be necessary to "float" the lip 
shock rather than smear it over several mesh points and thus destroy the 
detail of the flow in this region. 
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Concluding Remarks 

Solutions of the Navier-Stokes equations for the complete flow 
surrounding blunt bodies in a supersonic stream can be achieved using 
the numerical method described herein. Good comparisons with 
experimental data support this statement. The given coordinate system 
is especially useful in that it allows one to study the differences 
in the wake flow structure for many realistic body shapes. For 
example, the Viking Aeroshell exhibited a double recirculating flow 
pattern, while at similar freestream conditions, a sphere exhibited only 
a single pattern. However, it is evident that certain areas of the 
wake flow will require special treatment for increasing Reynolds 
number (i.e., Reynolds numbers on the order of 10,000). In particular, 
for blunt bodies experiencing a rapid expansion around a sharp corner 
(i.e. Viking Aeroshell) it appears that a lip shock will need to be 
treated as a discontinuity. The recompression shock may also need 
to be treated as a discontinuity and the mesh resolution in its 
vicinity to be increased in order to improve calculations in the far 
wake (i.e., 3 ^ x < lOO) . 
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Table 1 Calculated properties of Viking flowfields 


M 

Re 

Co 

deg 
sep ^ 

L 

Pw(o)/p^ 

Tw(o)/T„ 

Ptnax^Poo 

2 

1,000 

1.574 

80.6 

3.82 

0.737 

1.53 

1.090 

2 

5,000 

1.565 

88.2 

3.34 

0.901 

1.55 

1.152 

2 

30,000 

1.566 

93.6 

3.36 

0.943 

1.56 

1.179 

5 

5,000 

1.532 

72.0 

2.64 

1.34 

2.85 

2.21 

10 

1,000 

1.520 

25.2 

0.72 

3.08 

12.3 

4.88 

10 

5,000 

1.543 

55.1 

2.6 

3.63 

13.2 

6.15 
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Fig. 1. Analytic approximation to Viking Aerosheli with associated 
coordinate lines. 
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Fig. 3, Ori9nt3tion of vGlocity voctors in transformed coordinate system. 
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Fig. 4. Distribution of Cp and p/p^ over cylinder, Tq-300°K. 
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Fig. 5. Distribution of q over cylinder; M-5.5, Re * 1085, Tw/To“0.3, Tq-300°K 

M-5.73, Re-2050. Tw/To-0.68, To-300“ K 
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Fig. 6. Static pressure distribution on wake centerline 



Fig, 7. Velocity vectors over cylinder, iVi^-5o64, Rea,» 8350o 




Fig. 8. Velocity vectors over sphere, Mg^*5.64, Re^* 8350. 







Fig. 11. Velocity vectors over Viking Aeroshell, M^«2, Rp -5000. 
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Fig. 12. Pressure distribution along wake symmetry axis for Viking Aerosheil 
/VU-2, Re^-5000. 
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